*------------------------------------------------------------------------------*
					*** FIGURE 1. POLITBAROMETER ***
*------------------------------------------------------------------------------*


***Prepare pollution and weather data

**KW


clear
set obs 33
gen year = 1988 + _n
expand 12
bys year: gen month = _n
expand 31
bys year month: gen day = _n

gen date = mdy(m,d,y)
format date %td

keep if date!=.

gen dow=dow(date)
recode dow (0=7)

sort dow date
gen double week_id = _n if dow==1


sort date

forval dow = 2/7 {
	replace week_id = week_id[_n-`dow'+1] if dow==`dow'
} // dow

bys year dow (date): gen kw = _n if dow==4

bys week_id: egen kw2 = total(kw), miss
bys week_id: egen date_mon = total(date*(dow==1)), miss

format date_mon %td

drop kw 
rename kw2 kw


bys kw year dow : keep if _N==1

sort date

qui compress
save "${tempdata}kw.dta", replace

*****************************************

**Pollution

use "${rawdata}pollution/pollution_UBA_d_19902019.dta", clear

keep if year >= 2000

gen date = mdy(m,d,y)
format date %td

merge m:1 date using "${tempdata}kw.dta", keep(master match) nogen keepusing(kw week_id)
 
gen bula = substr(station,3,2)
 
gen bl 		= 1 if bula=="SH" 
replace bl 	= 2 if bula=="HH" 
replace bl 	= 3 if bula=="NI" 
replace bl 	= 4 if bula=="HB" 
replace bl 	= 5 if bula=="NW" 
replace bl 	= 6 if bula=="HE" 
replace bl 	= 7 if bula=="RP" 
replace bl 	= 8 if bula=="BW" 
replace bl 	= 9 if bula=="BY" 
replace bl 	= 10 if bula=="SL" 
replace bl 	= 11 if bula=="BE" 
replace bl 	= 12 if bula=="BB" 
replace bl 	= 13 if bula=="MV" 
replace bl 	= 14 if bula=="SN" 
replace bl 	= 15 if bula=="ST" 
replace bl 	= 16 if bula=="TH" 

drop if bl==.

collapse (mean) PM10, by(bl year kw week_id) 

bys bl year kw: drop if _N>1
bys bl week_id: drop if _N>1

qui compress
save "${tempdata}PM.dta", replace

*****************************************

**Weather

use "${rawdata}pollution/weather_d_19802019.dta", clear


gen year = year(date)

keep if year >= 2000

merge m:1 station 	using "${rawdata}pollution/meta_station_DWD.dta", keep(master match) nogen keepusing(bula)
merge m:1 date 		using "${tempdata}kw.dta", keep(master match) nogen keepusing(kw week_id)
 
gen bl 		= 1 if bula=="Schleswig-Holstein" 
replace bl 	= 2 if bula=="Hamburg" 
replace bl 	= 3 if bula=="Niedersachsen" 
replace bl 	= 4 if bula=="Bremen" 
replace bl 	= 5 if bula=="Nordrhein-Westfalen" 
replace bl 	= 6 if bula=="Hessen" 
replace bl 	= 7 if bula=="Rheinland-Pfalz" 
replace bl 	= 8 if bula=="Baden-Württemberg" 
replace bl 	= 9 if bula=="Bayern" 
replace bl 	= 10 if bula=="Saarland" 
replace bl 	= 11 if bula=="Berlin" 
replace bl 	= 12 if bula=="Brandenburg" 
replace bl 	= 13 if bula=="Mecklenburg-Vorpommern" 
replace bl 	= 14 if bula=="Sachsen" 
replace bl 	= 15 if bula=="Sachsen-Anhalt" 
replace bl 	= 16 if bula=="Thüringen" 

drop if bl==.

collapse (mean) temp wind_speed precip_amount, by(bl year kw week_id) 

bys bl year kw: drop if _N>1
bys bl week_id: drop if _N>1



	qui compress
	save "${tempdata}weather.dta", replace




*------------------------------------------------------------------------------*

***Politbarometer data

*This we cannot share
use "${rawdata}politbarometer/ZA2391_v12-0-0.dta", clear

rename v3  month
rename v4  year
rename v80 kw
rename v75 bl


keep if year >= 2003


rename v15 appr_gov
rename v16 appr_opp
rename v19 appr_elite
rename v21 interest_pol
rename v22 spectrum_lr
recode v25 (1 2=1) (3 4 5 = 0) (0 9 = .), gen(econ_good)
recode v25 (4 5=1) (1 2 3 = 0) (0 9 = .), gen(econ_bad)
recode v25 (3=1) (1 2 4 5 = 0) (0 9 = .), gen(econ_neutral)
recode v27 (1 2=1) (3 4 5 = 0) (0 9 = .), gen(ownecon_good)
recode v27 (4 5=1) (1 2 3 = 0) (0 9 = .), gen(ownecon_bad)
recode v27 (3 =1) (1 2 4 5= 0) (0 9 = .), gen(ownecon_neutral)
rename v54 gender
rename v56 age
rename v57 marst
rename v60 educ
rename v64 emplst
recode v65 (5/8=2) (9/16=3) (17/20=4) (21=1) (24/25=5) (0 99 = .), gen(occst)
rename v66 hhsize
recode v71 (0 9 = .) (4 = 0), gen(jobatrisk)
recode v72 (1/9 11 = 1 "yes") (10 98 =0 "no") (99=.), gen(leantoparty)
rename v77 citysize
rename V81 weight

*Sonntagsfrage
gen CDU 		= v6==1 	if v6<96
gen SPD 		= v6==2 	if v6<96
gen FDP 		= v6==3 	if v6<96
gen GRUNE 		= v6==4 	if v6<96
gen DIELINKE 	= v6==6 	if v6<96
gen AfD 		= v6==49 	if v6<96 

gen dow = 5

merge m:1 kw year dow using "${tempdata}kw.dta", keep(master match) nogen keepusing(date_mon)

*Generate Incumbent
gen Incumbent = .
label variable Incumbent "Incumbent"

	//BW
	replace Incumbent = ( SPD + FDP ) 		if Incumbent==. & date_mon<= td(05oct1980)
	replace Incumbent = ( CDU + FDP ) 		if Incumbent==. & date_mon<= td(06mar1983)
	replace Incumbent = ( CDU + FDP ) 		if Incumbent==. & date_mon<= td(25jan1987)
	replace Incumbent = ( CDU + FDP ) 		if Incumbent==. & date_mon<= td(02dec1990)
	replace Incumbent = ( CDU + FDP ) 		if Incumbent==. & date_mon<= td(16oct1994)
	replace Incumbent = ( CDU + FDP ) 		if Incumbent==. & date_mon<= td(27sep1998)
	replace Incumbent = ( SPD + GRUNE ) 	if Incumbent==. & date_mon<= td(22sep2002)
	replace Incumbent = ( SPD + GRUNE ) 	if Incumbent==. & date_mon<= td(18sep2005)
	replace Incumbent = ( CDU + SPD) 		if Incumbent==. & date_mon<= td(27sep2009)
	replace Incumbent = ( CDU + FDP ) 		if Incumbent==. & date_mon<= td(22sep2013)
	replace Incumbent = ( CDU + SPD ) 		if Incumbent==. & date_mon<= td(24sep2017)
	replace Incumbent = ( CDU + SPD ) 		if Incumbent==. & date_mon<= td(24sep2021)
	
*Generate Established opposition
gen Est_opp = .
label variable Est_opp "Established opposition"

	//BW
	replace Est_opp = ( CDU + GRUNE + DIELINKE ) 					if Est_opp==. & date_mon<= td(05oct1980)
	replace Est_opp = ( SPD + GRUNE + DIELINKE ) 					if Est_opp==. & date_mon<= td(06mar1983)
	replace Est_opp = ( SPD + GRUNE + DIELINKE) 					if Est_opp==. & date_mon<= td(25jan1987)
	replace Est_opp = ( SPD + GRUNE + DIELINKE) 					if Est_opp==. & date_mon<= td(02dec1990)
	replace Est_opp = ( SPD + GRUNE + DIELINKE) 					if Est_opp==. & date_mon<= td(16oct1994)
	replace Est_opp = ( SPD + GRUNE + DIELINKE) 					if Est_opp==. & date_mon<= td(27sep1998)
	replace Est_opp = ( CDU + FDP + DIELINKE ) 						if Est_opp==. & date_mon<= td(22sep2002)
	replace Est_opp = ( CDU + FDP + DIELINKE ) 						if Est_opp==. & date_mon<= td(18sep2005)
	replace Est_opp = ( FDP + GRUNE + DIELINKE ) 					if Est_opp==. & date_mon<= td(27sep2009)
	replace Est_opp = ( SPD + GRUNE + DIELINKE) 					if Est_opp==. & date_mon<= td(22sep2013)
	replace Est_opp = ( FDP + GRUNE + DIELINKE ) 					if Est_opp==. & date_mon<= td(24sep2017)
	replace Est_opp = ( FDP + GRUNE + DIELINKE + AfD ) 				if Est_opp==. & date_mon<= td(26sep2021)

* Other parties	
gen Other = 0 if v6<96
replace Other = 1 if Other==0 & Incumbent==0 & Est_opp==0 

* Drop Sonntagsfrage for interviews AFTER election day in election years
gen post_election = 0 

replace post_election = 1 if year==1980 & date_mon > td(05oct1980)
replace post_election = 1 if year==1983 & date_mon > td(06mar1983)
replace post_election = 1 if year==1987 & date_mon > td(25jan1987) 
replace post_election = 1 if year==1990 & date_mon > td(02dec1990) 
replace post_election = 1 if year==1994 & date_mon > td(16oct1994) 
replace post_election = 1 if year==1998 & date_mon > td(27sep1998)
replace post_election = 1 if year==2002 & date_mon > td(22sep2002)
replace post_election = 1 if year==2005 & date_mon > td(18sep2005)
replace post_election = 1 if year==2009 & date_mon > td(27sep2009)
replace post_election = 1 if year==2013 & date_mon > td(22sep2013)
replace post_election = 1 if year==2017 & date_mon > td(24sep2017)

replace Incumbent 	= . if post_election==1
replace Est_opp 	= . if post_election==1
replace Other 		= . if post_election==1

rename Incumbent 	sonntag_incumb
rename Est_opp 		sonntag_oppos
rename Other 		sonntag_other

drop v* V*

recode appr_gov appr_opp emplst spectrum_lr (0 99 = .)
recode appr_elite (0 9 = .) (2 3 = 0)

recode appr_gov (7/11 = 1) (1/6 = 0), gen(appr_gov_pos)
recode appr_gov (6/11 = 0) (1/5 = 1), gen(appr_gov_neg)
recode appr_gov (6 = 1) (1/5 7/11= 0), gen(appr_gov_neutral)

recode appr_opp (7/11 = 1) (1/6 = 0), gen(appr_opp_pos)
recode appr_opp (6/11 = 0) (1/5 = 1), gen(appr_opp_neg)
recode appr_opp (6 = 1) (1/5 7/11= 0), gen(appr_opp_neutral)

gen appr_diff = appr_gov - appr_opp

gen appr_gov_gt_opp = appr_gov > appr_opp if appr_gov!=. & appr_opp!=.
gen appr_gov_lt_opp = appr_gov < appr_opp if appr_gov!=. & appr_opp!=.

gen female = gender==2

recode age hhsize bl (99=.)
recode spectrum_lr (99 0 =.) (1/5=1) (6/11=0), gen(leftofcent)
recode spectrum_lr (99 0 =.) (1/6=0) (7/11=1), gen(rightofcent)
recode age (1 2 3 = 1 "age<30") (4/8=2 "age 30-59") (9/11 = 3 "age 60+") (99=.), gen(age2)
recode marst (0 9 =.) (2 = 1)
recode citysize (9=.) (1 2 3 = 4) (5=6) (8=7)
recode educ (5=1) (6=7) (9=.)
recode interest_pol (0 8 =.)
recode interest_pol (1/2=1) (3/5=0), gen(interest_pol_high)
recode interest_pol (3=1) (1/2 4/5=0), gen(interest_pol_neutral)
recode interest_pol (1/3=0) (4/5=1), gen(interest_pol_low)

replace bl = bl-1 if bl>=12

merge m:1 bl year kw using "${tempdata}PM.dta", keep(master match) nogen 
merge m:1 bl year kw using "${tempdata}weather.dta", keep(master match) nogen 

*Lags and leads

gen PM10_gt20 = PM10_ug > 20 if PM10_ug!=.

rename PM10_ug 		PM10_orig
rename PM10_gt20 	PM10_gt20_orig

rename week_id week_id_orig

forval l = 1/2 {

	//Lags
	gen week_id = week_id_orig - `l'	
	merge m:1 bl week_id using "${tempdata}PM.dta", 	keep(master match) nogen keepusing(PM10_ug)		
	gen PM10_gt20_lag`l' = PM10_ug > 20 if PM10_ug!=.	
	rename PM10_ug PM10_lag`l'	
	drop week_id
	
	//Leads
	gen week_id = week_id_orig + `l'	
	merge m:1 bl week_id using "${tempdata}PM.dta", 	keep(master match) nogen keepusing(PM10_ug)	
	gen PM10_gt20_lead`l' = PM10_ug > 20 if PM10_ug!=.
	rename PM10_ug PM10_lead`l'	
	drop week_id	

}  // l

rename week_id_orig week_id


qui compress
save "${tempdata}estimation.dta", replace




*------------------------------------------------------------------------------*

***Regressions and graph


use "${tempdata}estimation.dta", clear


*****************************************

**Regressions

egen clust = group(bl year)

recode PM10_orig (min/15=0) (15/20=1) (20/25=2) (25/30=3) (30/max=4), gen(PM10_orig_dum)

xtile PM10_orig_q4 = PM10_orig, n(4)

gen age60 = age2==3 if age2<.

local outcomes  	"sonntag_* appr_*_neg appr_*_pos appr_*_neutral appr_elite ownecon* econ* leantoparty interest_pol_*"
local treatments 	"PM10_gt20_orig"
local fe 			"year#kw bl#month age female educ city marst emplst occst" //"
local controls 		"c.temp##c.temp c.wind_speed##c.wind_speed c.precip_amount##c.precip_amount"

postutil clear	
postfile results str25(outcome treatment regression) double(beta cil cih N_obs mean_outcome sd_outcome mean_PM sd_PM Fstat) using "${tempdata}results.dta" , replace
	

foreach out of varlist `outcomes'  {

	foreach PM of varlist `treatments' {		
	
		local l = subinstr("`PM'","PM10_gt20_","",1)
		local l = subinstr("`l'","PM10_","",1)
	
		foreach reg in OLS {
	
			if "`reg'"=="OLS" eststo `out': reghdfe `out' `PM'	`controls' [aw=weight] , a(`fe') vce(cluster clust) 
			
			qui sum `out' if e(sample)==1
			local mean_out 	= r(mean)
			local sd_out 	= r(sd)
			
			qui sum `PM' if e(sample)==1
			local mean_PM 	= r(mean)
			local sd_PM 	= r(sd)
			
			*Treatment coefficients
			cap di _b[`PM']
																	
			if _rc==0 {
				local b   	= _b[`PM'] 
				local cil   = _b[`PM'] - invt(e(df_r),0.975) * _se[`PM']
				local cih   = _b[`PM'] + invt(e(df_r),0.975) * _se[`PM']	
				local N 	= e(N)
				
				if "`reg'"=="IV" local F = e(widstat)
													
				if "`reg'"=="OLS" post results ("`out'") ("`PM'") ("`reg'") (`b') (`cil') (`cih') (`N') (`mean_out') (`sd_out') (`mean_PM') (`sd_PM') (.)
				if "`reg'"=="IV"  post results ("`out'") ("`PM'") ("`reg'") (`b') (`cil') (`cih') (`N') (`mean_out') (`sd_out') (`mean_PM') (`sd_PM') (`F')
				
			} //if
		
		} //reg
	
	} //PM
	
} //out
	
postclose results	

*****************************************

**Coef plots

local PM "PM10_gt20_orig"

coefplot ///
	sonntag_incumb, bylabel("Government parties") || ///
	sonntag_oppos, bylabel("Opposition parties") || /// 
	sonntag_other, bylabel("Other parties") ||  ///
	appr_gov_neg, bylabel("negative") ||  ///
	appr_gov_neutral, bylabel("neutral") ||  ///
	appr_gov_pos, bylabel("positive") ||  ///
	appr_opp_neg, bylabel("negative") ||  ///
	appr_opp_neutral, bylabel("neutral") ||  ///
	appr_opp_pos, bylabel("positive") ||  ///
	appr_elite, bylabel("positive") || ///
	, keep(`PM') xline(0, lc(cranberry)) xtitle("Effect of {bf:1}(PM10>20)") grid(none)  ///
	headings( ///
	1="{bf:Voting intention}" ///
	4="{bf:Approval of government}" ///
	7="{bf:Approval of opposition}"  ///
	10="{bf:Approval of elites}") /// 
	bycoefs	 scheme(lean2) graphregion(color(white)) color(navy) ciopts(lc(navy)) ms(oh) swapnames


graph export "${outputs}figure1.png", replace width(3000)


*------------------------------------------------------------------------------*

clear

exit
